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We present a unified framework for renormalization group methods, including Wilson's numerical 
renormalization group (NRG) and White's density-matrix renormalization group (DMRG), within 
the language of matrix product states. This allows improvements over Wilson's NRG for quantum 
impurity models, as we illustrate for the one-channel Kondo model. Moreover, we use a variational 
method for evaluating Green's functions. The proposed method is more flexible in its description 
of spectral properties at finite frequencies, opening the way to time-dependent, out-of-equilibrium 
impurity problems. It also substantially improves computational efficiency for one-channel impurity 
problems, suggesting potentially linear scaling of complexity for n-channel problems. 
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Wilson's numerical renormalization group (NRG) is a 
key method 0] for solving quantum impurity models such 
as the Kondo, Anderson or spin-boson models, in which 
a local degree of freedom, the "impurity" , is coupled to 
a continuous bath of excitations. These models are of 
high relevance in the description of magnetic impurities, 
of quantum dots, and problems of decoherence. NRG 
has been used with great success to calculate both ther- 
modynamic 0, and dynamical 0, i, i] properties. It 
is, however, of limited use in more complex situations: 
Computational cost grows exponentially for a coupling 
to multiple bands in the bath. In systems out of equilib- 
rium or with time-dependent external parameters, such 
as occur in the tuning of quantum dots, difficulties arise 
due to NRG's focus on low energy properties through its 
logarithmic discretization scheme which looses accuracy 
at high spectral frequencies. 

In the present Letter, we draw attention to the fact 
that states generated by the NRG have the structure of 
matrix product states (MPS) 0, Q on a 1-D geometry. 
This is a simple observation, which however has impor- 
tant conceptual and practical implications: 

(i) As White's density matrix renormalization group 
(DMRG) [| for treating quantum chain models is in its 
single-site version identical to variational MPS Q, NRG 
and DMRG are now seen to have the same formal basis 
of matrix product states, resolving a long-standing ques- 
tion about the connection between both methods, (ii) 
All NRG results can be improved upon systematically 
by variational optimization in the space of variational 
matrix product states (VMPS) of the same structure as 
those used by NRG. This does not lead to major improve- 
ments at us — where NRG works very well, but leads to 
the inclusion of feedback from low- to high-energy states, 
also allowing the relaxation of the logarithmic bath dis- 



cretization of NRG: spectra away from to = can be 
described more accurately and with higher resolution, 
(iii) Recent algorithmic advances using VMPS Q An par- 
ticular those treating time-dependent problems (ToL ll l| , 
can now be exploited to tackle quantum impurity mod- 
els involving time-dependence or nonequilibrium; this in- 
cludes applications to the description of driven qubits 
coupled to decohering baths, as relevant in the field of 
quantum computation, (iv) The VMPS algorithm allows 
ground state properties of quantum impurity models to 
be treated more efficiently than NRG: the same accu- 
racy is reached in much smaller ansatz spaces (roughly, of 
square-root size). Moreover, our results suggest that for 
many (if not all) rt-channcl impurity problems it should 
be feasible to use an unfolded geometry, for which the 
complexity will only grow linearly with n. 

The present Letter provides a "proof of principle" for 
the VMPS approach to quantum impurity models by ap- 
plying it to the one-channel Kondo model. We repro- 
duce the NRG flow of the finite size spectrum [2j, and 
introduce a VMPS approach for calculating Green's func- 
tions, as we illustrate for the impurity spectral function 
[3[, which yields a significant improvement over existing 



alternative techniques [12], [13|, [14|, [15| . Our results illus 



trate in which sense the VMPS approach is numerically 
more efficient than the NRG. 

NRG generates matrix product states: — To be specific, 
we consider Wilson's treatment of the Kondo model, de- 
scribing a local spin coupled to a fermionic bath. To 
achieve a separation of energy scales, the bath excita- 
tions are represented by a set of logarithmically spaced, 
discrete energies um = A~™, where A > 1 is a "discretiza- 
tion parameter" By tridiagonalization, the model 
is then mapped onto the form of a semi-infinite chain 
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H = limjv_ 00 H JV where [l| 

N-l 

H N = -2Js -S+Y,tn (c^c n+1 ,^ + (!) 

71=1 

with S = ^CqCtco and creation (annihilation) operators 
(c) , respectively. 7^ describes an impurity spin s cou- 
pled to the first site of a chain of length N of Fermions 
with spin /i and exponentially decreasing hopping ma- 
trix elements along the chain (£„ ~ A~™/ 2 ). TL N lives on 
a Hilbert space spanned by the set of did N basis states 
{\io, ii, «2, ■ ■ • *jv)}, where io labels the di possible impu- 
rity states and i n (for n — 1, . . . , N) the d possible states 
of site n (for the Kondo model, io = {T,4} and for all 
other sites i n = {0, f, J,, f J,}, i. c. d j = 2 and d = 4). 

To diagonalize the model, NRG starts with a chain 
of length (n — 1), chosen sufficiently small that 7i™ -1 
can be diagonalized exactly, yielding a set of eigenstates 
|^>™ -1 ). One continues with the subsequent iterative pre- 
scription: project H. n ~ 1 onto the subspace spanned by 
its lowest D eigenstates, where D < did 71-1 is a control 
parameter (typically between 500 and 2000); add site n 
to the chain and diagonalize Ti n in the enlarged (Dd)- 
dimcnsional Hilbert space, writing the eigenstates as 

if ^Eirt^S 1 , (2) 

in—I a—1 

where the coefficients have been arranged in a matrix 
Paj^ with matrix indices a, /3, labelled by the site index 
n and state index i n ; rescale the eigenenergies by a factor 
A 1 / 2 ; and repeat, until the eigenspectrum converges, typ- 
ically for chain lengths N of order 40 to 60. At each step 
of the iteration, the eigenstates of TL N can thus be writ- 
ten [by repeated use of Eq. in the form of a so-called 
matrix product state, 

IO = P% ] P&l 1 P^---P&lJio,ii,..-,iN) (3) 

(summation over repeated indices implied). The ground 
state is then the lowest eigenstate of the effective Hamil- 
tonian Ti^p = i' l Pa\'^- N \i'p)^ 1 - e - the projection of the 
original Ti on the subspace of MPS of the form ||3J). 

VMPS optimization: — Let us now be more ambitious, 
and aim to find the best possible description of the ground 
state within the space of all MPSs of the form ([3]), using 
the matrix elements of the matrices {p[ n l} with P<- n ' = 
{P''"'} as variational parameters to minimize the energy. 
Using a Lagrange multiplier to ensure normalization, we 
thus study the following optimization problem: 

min [(tp N \H N \ij N )-\(iP N \^ N )] . (4) 
| 1 /- N )e{MPS} L 

This cost function is multiquadratic in the dj + d(N — 1) 
matrices {P^™'} with a multiquadratic constraint. Such 



problems can be solved efficiently using an iterative 
method in which one fixes all but one (let's say the n'th) 
of the matrices {Pl n l} at each step; the optimal P^ min- 
imizing the cost function given the fixed values of the 
other matrices can then be found by solving an eigen- 
value problem ||. With P^ optimized, one proceeds 
the same way with pl rl + 1 ] and so on. When all matrices 
have been optimized locally, one sweeps back again, and 
so forth. By construction, the method is guaranteed to 
converge as the energy goes down at every step of the it- 
eration, having the ground state energy as a global lower 
bound. Given the rather monotonic hopping amplitudes, 
we did not encounter problems with local minima. 

In contrast, NRG constructs the ground state in a sin- 
gle one-way sweep along the chain: each is thus 
calculated only once, without allowing for possible feed- 
back of P's calculated later. Yet viewed in the above 
context, the ground state energy can be lowered further 
by MPS optimization sweeps. This accounts for the feed- 
back of information from low to high energy scales. This 
feedback may be small in practice, but it is not strictly 
zero, and its importance increases as the logarithmic dis- 
cretization is refined by taking A — > 1. Note that the 
computational complexity of both VMPS optimization 
and NRG scales as NdD 3 @, 0| , and symmetries can be 
exploited (with similar effort) in both approaches. The 
inclusion of feedback leads to a better description of spec- 
tral features at high frequencies, which are of importance 
in out-of-equilibrium and time-dependent impurity prob- 
lems. Moreover, it also allows to relax the logarithmic 
discretization scheme, further improving the description 
of structures at high frequency as illustrated below. 

Energy level flow: — The result of a converged set of 
optimization sweeps is a VMPS ground state IV'o'Lpf the 
form ©; exploiting a gauge degree of freedom [8(, the 
P's occurring therein can always be chosen such that all 
vectors = [pi 10 ! . . . pl l "l] |j 0i are orthonor- 

mal. The effective Hamiltonian at chain length n, the 
central object in NRG, is then H n af} = {^\A n / 2 H n \'ip^) . 
Its eigenspectrum can be monitored as n increases, re- 
sulting in an energy level flow along the chain. 

Green's functions: — Similar techniques also allow 
Green's functions to be calculated variationally [l5| . The 
typical Green's functions of interest are of the form 
G°(ui) = (it>p\c\\) where \x), commonly called a correc- 
tion vector [17| . is defined by 

|X> = — c + |^o> , (5) 

u> — ri + in 

with |?/?o) the ground state of the system, e. g. calcu- 
lated using the VMPS approach and thus represented 
as MPS. The spectral density is then given by A(u) — 
— lim^—vo jSm [G^ (w)) . The (unnormalized) state |x) 
may be calculated variationally within the set of MPS 



3 



by optimizing the weighted norm 
1 
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W=(H-u) 2 +v 2 



where \\0\\w = (€\ W \0> and weight W > such that 
it yields a quadratic equation. Minimizing M efficiently 
by optimizing one P at a time leads to two independent 
optimizations over 3?e|x) and 9ta|x), respectively. Both 
involve only multilinear terms such that each iteration 
step requires to solve a sparse linear set of equations [l(| • 
The minimization of Af in Eq. ([6|), however, can be- 
come increasingly ill-conditioned as r\ — > 19], with con- 
ditioning deteriorating quadratically in r\. If one directly 
solves 6/6(PW I [(xKft-c- i v )\ x ) _ (x\J\ip)] = by 
a non-hermitian equation solver such as the biconjugate 
gradient method, conditioning deteriorates only linearly. 

Application to Kondo model: — Let us now illustrate 
above strategies by applying them to the Kondo model. 
Since the Hamiltonian in Eq. {1} couples f and [ band 
electrons only via the impurity spin, it is possible (see 
also H EH) to "unfold" the semi-infinite Wilson chain 
into an infinite one, with | band states to the left of 
the impurity and J, states to the right, and hopping am- 
plitudes decreasing in both directions as A _ l"l/ 2 . Since 
the left and right end regions of the chain, which de- 
scribe the model's low-energy properties, are far apart 
and hence interact only weakly with each other (ana- 
lyzed quantitatively in terms of mutual information in 
Fig. [TJd), the effective Hamiltonian for these low ener- 
gies will be of the form Tif <£> t L + 1 T ® Tif. Due to 
the symmetry of the Kondo coupling, 7i| ff and 7iJ ff have 
the same eigenspectrum for n » 1, such that the fixed 
point spectrum is already well-reflected by analyzing ei- 
ther one, as illustrated in Fig.QJa). Note that for a direct 
comparison with NRG, the spin chains can be recombined 
within VMPS 
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The resulting standard energy flow 
diagram presented in panel (a) for VMPS and NRG, re- 
spectively, show excellent agreement for low energies for 
all n including the fixed point spectrum. 

The dimensions of the effective Hilbert spaces needed 
for VMPS and NRG to capture the low energy properties 
(here energy resolution better than Tk) are roughly re- 
lated by Dmps ~ V -Dnrg [111)' implying significant com- 



putational gain with VMPS, as calculation time scales 
as D 3 for both. Indeed, Fig. [He) shows that VMPS has 
three orders of magnitude of better precision for the same 
D. More generally, if the impurity couples to n electronic 
bands (channels) , the Wilson chain may be unfolded into 
a star-like structure of 2n branches, with Z?mps ~ D 



l/2n 
NRG' 



This implies that for maintaining a desired precision in 
going from 1 to n channels, -Dmps will stay roughly con- 
stant, and calculation time for all sites other than the 
impurity will scale merely linearly with the number of 
channels. Whether the chains can be unfolded in prac- 
tice can easily be established by checking whether or not 




f-. 


VMPS (D=32) 




— odd 




\, • even 


(c) 





10 20 30 
Wilson shell n 



16 32 64 128 256 512 
dimension D 



Figure 1: (Color online) Comparison of VMPS and NRG 
data for logarithmic discretization, (a) Energy level flow of 
the Kondo model as a function of site index n obtained from 
Ti^f of a variationally optimized MPS with Dmps = 32 (light 
red), of the corresponding recombined spin chains (red) jla ]. 
and from NRG using Dnrg = 32 2 states (dashed black), 
(b) Correlation along the Wilson chain between spin-up and 
spin-down at site n in terms of mutual information Jm(w) = 
S(n^) + S(ni) — S(n^,ni). Here S is the entropy of the re- 
duced density matrix of the groundstate with respect to the 
indicated subspace [l|| (solid for even, dashed for odd sites 
n). The Wilson shell corresponding to Tk is indicated by the 
vertical dashed line, (c) Bond entropy S along the unfolded 
Wilson chain where S is the usual von Neumann entropy of 
the VMPS reduced density matrix when going from site n to 
n + 1. (d) Comparison of T-matrix (SJmT^, see also Fig. [2} 
for B = between VMPS and NRG, including deconvoluted 
VMPS data [ijj. Inset shows zoom into peak at uj = 0. The 
significantly smaller A = 1.2 applicable for VMPS (discretiza- 
tion intervalls are indicated by vertical lines) shows clearly im- 
proved agreement with the Friedel sum rule T (0) 7r 2 /2 = 1. 
(e) Comparison of ground state energy of the Kondo Hamil- 
tonian {TJ for fixed chain length relative to the extrapolated 
energy for D — > oo for VMPS and NRG as function of the 
dimension D of states kept. 



the correlation between them, characterized e.g. in terms 
of mutual information, decays rapidly with increasing n 
(cf. Fig. [TJd and caption) . 

Adaptive discretization: — Through its variational 
character, VMPS does not rely on logarithmic discretiza- 
tion crucial for NRG. The potential of greatly enhanced 
energy resolution using VMPS is already indicated by 
the A = 1.2 data in Fig. [[Yd). It is illustrated to full 
extent in Fig. [2l showing the splitting of the Kondo peak 
in the presence of a strong magnetic field calculated us- 
ing VMPS (bare: dots, deconvoluted: red solid), stan- 
dard NRG (blue dashed) and perturbatively [2l[ (black). 
By using a linear (logarithmic) discretization scheme for 
|u;| <B (\uj \ > B), respectively, VMPS yields well-resolved 
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Figure 2: (Color online) Impurity spectral function for the 
Kondo model -N F QmT^(uj) = J 2 {{0\ |C M >)„ for B > Tk , 
where C M = S ■ cr^ic , and Nf is the density of states at 
the Fermi energy, calculated with VMPS (dots: raw data, red 
solid: deconvoluted) , NRG (dashed) and perturbative (black 
solid) [2l| . According to [2l|. the peak of the perturbative 
result should be shifted in uj by — B /2\og(B /Tk) (arrow). 
NRG and VMPS discretization intervals are indicatd by blue 
shaded areas and gray vertical lines, respectively. The in- 
set shows the hopping amplitudes corresponding to standard 
(blue dashed) and adapted (red solid) discretization schemes. 
The required Lorentzian broadening rj of the VMPS data 
smears out sharper features. Deconvolution (targeting with 
adaptive spline) together with subsequent Gaussian broaden- 
ing was applied to obtain the red solid line fl9( ] . 



sharp spectral features at finite frequencies. This resolu- 
tion is out of reach for NRG, whose discretization inter- 
vals (blue shaded intervals), even for comparatively small 
choice of A = 1.7, are much broader than the spectral fea- 
tures of interest. The line shape of our deconvoluted data 
(red solid line) agrees well with the analytic RG calcula- 
tion [2l[ (black solid line), perturbative in 1/ \og(B /Tk)- 
The peak positions agree well also after a shift in u> by 
—B/2 log(_B/Tx) of the perturbative result suggested by 
PH is taken into account. 

Outlook: — Let us finish by pointing out that the MPS 
approach can readily be extended to the case of finite 
temperatures by using matrix product density operators 
fioj ] instead of MPS, and to time-dependent problems 
(such as TL = Ti{t) or non-equilibrium initial conditions), 
by usin g th e recently developed adaptive time-dependent 
DMRG [lH and MPS analogues thereof [UJ. Preliminary 
work in this direction was very encouraging and will be 
published in the near future. 

In conclusion, the MPS approach provides a natural 
language for simulating quantum impurity models. The 
underlying reason is that these models, when formulated 
on the Wilson chain, feature only nearest-neighbor in- 
teractions. Their low-energy states are thus determined 



mainly by their nearest-neighbor reduced density matri- 
ces, for which very good approximations can be obtained 
by suitably optimizing the set of matrices constituing a 
MPS [2(J • We also showed how these could be used for a 
direct variational evaluation of Green's functions. 
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APPENDIX - DECONVOLUTION OF 
SPECTRAL DATA 

DMRG obtains spectral data from a discretized model 
Hamiltonian. In order for the spectral data to be smooth, 
an intrinsic frequency dependent Lorentzian broadening 
r\ is applied during the calculation of the correction vector 
\x)k at frequency cu k (cf. Eq. [5]), 

S Vk (oj - oj k ) = — - — j • (7) 

7T (oj - UJ k ) Z + 1]* 

Since the original model has a continuous spectrum, the 
broadening r\ k should be chosen of the order or larger 
than the artificial coarse grained discretization intervals 
8 U . Larger rj of course improves numerical convergence. 
However, since Lorentzian broadening produces longer 
tails than for example Gaussian broadening, this makes it 
more susceptible to pronounced spectral features closeby. 
Our general strategy for more efficient numerical treat- 
ment was then as follows, (i) Choose somewhat larger r\ 
(j\ ~ 28 u) throughout the calculation, (ii) Deconvolve the 
raw data to such an extent that the underlying discrete 
structure already becomes visible again, (iii) followed by 
a Gaussian smoothening procedure which then acts more 
locally. Let us describe step (ii) in more detail. 

Broadening, by construction, looses information. 
Hence trying to obtain the original data from the 
broadened data via deconvolution is intrinsically ill- 
conditioned. In literature there are several ways of deal- 
ing with this problem, most prominently maximum en- 
tropy algorithms (see [B[ and reference below). Our ap- 
proach is targeting the actual spectrual function using 
the knowledge about the Lorentzian broadening used 
during the VMPS calculation, combined with adaptive 
spline. Given the data A (oj) obtained through VMPS, 
let us propose the existence of some smooth but a priori 
unknown target curve A (uj), which when broadened the 
same way as the VMPS data using exactly the same rj k 
via a Lorenzian broadening kernel 

oo 

A k = A(uj k )= J dJ A(uj')8 Vk (uj' -uj k ), (8) 

— oo 

reproduces the original data A(u>). Direct inversion of 
above equation as it is is ill-conditioned, as already men- 
tioned, and not useful in practice. 

Let us assume the unknown target curve A (oj) is 
smooth and parametrized by piecewise polynomials. 
Given the data points uj k with k = 1 . . . N , the inter- 
vals in between these values will be approximated in the 
spirit of adaptive spline functions by 3 rd order polynomi- 
als (k = 1 ... N - 1) 

( a k + b k (oj - uj k ) + c k (oj - uj k ) 2 + 
fk M = < d k (oj - oj k ) 3 for oj S [uj k , OJk+i] (9) 
I otherwise. 



Since spectral functions decay as 1/lu 2 for large oj, for 
our purposes the ends are extrapolated assymptotically 
to infinity, allowing both 1/oj and 1/uj 2 polynomials 

/o ^ 55 { otherwise 

fs M - ( + ^ "i ^ (10) 
J y ' \ otherwise. v ; 

In total, this results in 4 (N - 1) + 2 x 2 = AN pa- 
rameters, with the target function paramatrized piece- 
wise as A (oj) = f (oj) = J2 k =o f k ( w ) ■ I n cases where 
one has not approached the assymptotic limit yet, the 
ends may simply be modelled also by Eq. ([9|), taking 
Co = da = c/v = djsi — 0. Moreover, if information about 
the gradient f'(oj) is known, it can be built in straight- 
forwardly in the present scheme by replacing b k . 

The parameters for the piecewise parametrization are 
solved for by requiring the following set of conditions 

(i) The function / should be continuous and smooth 
by requiring that /, /' and /" are continuous (3N 
equations) . 

(ii) The function /, when broadened as in Eq. ([8]), should 

reproduce the VMPS data A k 

At^Y. I ^' 1* (-') , ; W " 2 ( n ) 

fe '=o / (w'-Wfc) +V k 

A k ~ A c k = Pk r k (12) 

where r k = f k (oj k ) - f^_i(<^k) and uj = -oo, 
lon+1 = +oo (N equations). 

In the spirit of adaptive spline, the third derivative of the 
piecewise polynomials is no longer required to be contin- 
uous. Its jump r k is set proportional to the change in 
A k — A c k introducing the additional prespecified parame- 
ter set pk, kept small for our purposes (note that enforc- 
ing the strict equality A c k = A k by setting p k = would 
result in an ill-conditioned problem). 

If intervall spacings specified by oj k are nonuniform, 
the p k have to be adapted accordingly. For this paper 
we used p k — p ■ (oJk+i — oj k ) a with p on the order of 
10~ 6 and a ~ 1. With p k fixed, Eqs. CP) and (jUJ) de- 
termine all spline parameters uniquely in terms of the 
original VMPS data A k . The integrals emerging out of 
Eq. (fTTj) can all be evaluated analytically. The final in- 
version of Eq. (fTTj) to obtain the parameters for f(oj) 
is well-behaved for small but finite p, small enough to 
clearly sharpen spectral features. 

Further reading 

W.H. Press, S.A. Teukolsky, W.T. Vetterling, B.P. Flan- 
nery, Numerical Recipies in C, 2nd ed., Cambridge Uni- 
versity Press, Cambridge (1993). 



